﻿/*
 * Coordinate
 * 
 * initial implementation : 12.12.2008
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 */
using System;
using System.Collections;
using System.Collections.Generic;
using Petra.Gis.Geometry;
using Petra.Gis.Geometry.Support;


namespace Petra.Gis.Geometry.Utils
{
  public static class SurfaceMath
  {



    public static bool IsCCW(LineRing geom)
    {
      ICoordinate hip; //, p;
      ICoordinate prev, next;
      int hii, i;
      int nPts = geom.coordList.Count;

      // check that this is a valid ring - if not, simply return a dummy value
      if (nPts < 4) return false;

      // algorithm to check if a Ring is stored in CCW order
      // find highest point
      //hip = geom.coordList[0];
      hii = 0;
      i = 0;

      double hiY = geom.coordList[0].Y;
      foreach (ICoordinate c in geom.coordList)
      {
        if (c.Y > hiY)
        {
          hiY = c.Y;
          hii = i;
        }
        i++;
      }

      hip = geom.coordList[hii];


      //for (i = 1; i < nPts; i++)
      //{
      //  p = geom.coordList[i];
      //  if (p.Y > hip.Y)
      //  {
      //    hip = p;
      //    hii = i;
      //  }
      //}

      // find points on either side of highest
      int iPrev = hii - 1;
      if (iPrev < 0) iPrev = nPts - 2;
      int iNext = hii + 1;
      if (iNext >= nPts) iNext = 1;
      prev = geom.coordList[iPrev];
      next = geom.coordList[iNext];

      // translate so that hip is at the origin.
      // This will not affect the area calculation, and will avoid
      // finite-accuracy errors (i.e very small vectors with very large coordinates)
      // This also simplifies the discriminant calculation.
      double prev2x = prev.X - hip.X;
      double prev2y = prev.Y - hip.Y;
      double next2x = next.X - hip.X;
      double next2y = next.Y - hip.Y;

      // compute cross-product of vectors hip->next and hip->prev
      // (e.g. area of parallelogram they enclose)
      double disc = next2x * prev2y - next2y * prev2x;
      // If disc is exactly 0, lines are collinear.  There are two possible cases:
      //	(1) the lines lie along the x axis in opposite directions
      //	(2) the line lie on top of one another
      //  (2) should never happen, so we're going to ignore it!
      //	(Might want to assert this)
      //  (1) is handled by checking if next is left of prev ==> CCW

      if (disc == 0.0)
      {
        // poly is CCW if prev x is right of next x
        return (prev.X > next.X);
      }
      else
      {
        // if area is positive, points are ordered CCW
        return (disc > 0.0);
      }

    }





    /// <summary>
    /// Returns the area of the LinearRing
    /// </summary>
    public static double GetArea(LineRing geom)
    {

        if (geom.coordList.Count < 3)
          return 0;

        double sum = 0;
        double ax = geom.coordList.startPoint.X;
        double ay = geom.coordList.startPoint.Y;

        double bx, by, cx, cy;

        IEnumerator<ICoordinate> pe = geom.coordList.GetEnumerator();
        pe.MoveNext();

        for (int i = 1; i < geom.coordList.Count - 1; i++)
        {
          //bx = geom.coordList[i].X;
          //by = geom.coordList[i].Y;
          //cx = geom.coordList[i + 1].X;
          //cy = geom.coordList[i + 1].Y;

          bx = pe.Current.X;
          by = pe.Current.Y;
          pe.MoveNext();
          cx = pe.Current.X;
          cy = pe.Current.Y;

          sum += ax * by - ay * bx +
            ay * cx - ax * cy +
            bx * cy - cx * by;
        }


        return Math.Abs(-sum / 2);
    }





    /// <summary>
    /// Returns true of the Point 'p' is within the instance of this ring
    /// </summary>
    /// <param name="p"></param>
    /// <returns></returns>
    public static bool IsPointWithin(LineRing geom, ICoordinate p)
    {
      bool c = false;
      for (int i = 0; i < geom.coordList.Count; i++)
        for (int j = i + 1; j < geom.coordList.Count - 1; j++)
        {
          if ((((geom.coordList[i].Y <= p.Y) && (p.Y < geom.coordList[j].Y)) ||
             ((geom.coordList[j].Y <= p.Y) && (p.Y < geom.coordList[i].Y))) &&
            (p.X < (geom.coordList[j].X - geom.coordList[i].X) * (p.Y - geom.coordList[i].Y) / (geom.coordList[j].Y - geom.coordList[i].Y) + geom.coordList[i].X))
            c = !c;
        }
      return c;
    }

  }
}
